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Abstract. 

Graphene placed in a magnetic field possesses an extremely high mid/far-infrarcd 
optical nonlinearity originating from its unusual band structure and selection rules for 
the optical transitions near the Dirac point. Here we study the linear and nonlinear 
optical response of graphene in strong magnetic and optical fields using quantum- 
mechanical density-matrix formalism. We calculate the power of coherent terahertz 
radiation generated as a result of four-wave mixing in graphene. We show that even 
one monolayer of graphene gives rise to appreciable nonlinear frequency conversion 
efficiency and Raman gain for modest intensities of incident infrared radiation. 
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1. Introduction and Background 

Graphene has unique electronic and optical properties stemming from linear, massless 
dispersion of electrons near the Dirac point and the chiral character of electron 
states[U |2]. Magnetooptical properties of graphene and thin graphite layers are 
particularly interesting, showing multiple absorption peaks and unique selection rules 
for transitions between Landau levels [31 HJ El E]- Recent progress in growing high- 
quality epitaxial graphene and graphite with high room-temperature mobility and strong 
magnetooptical response attracted a lot of interest and showed the promise of new 
applications in the infrared optics and photonics [3, El E] ■ The time is ripe to explore 
the nonlinear optical properties of a magnetized graphene and their applications. We 
have recently shown that graphene placed in a magnetic field possesses perhaps the 
highest infrared optical nonlinearity among known materials [5]. Here we present 
detailed derivation of the linear and nonlinear response of a magnetized graphene 
based on a rigorous density-matrix formalism. We apply this approach to calculate 
the terahertz radiation power generated by third-order nonlinear optical processes: the 
four-wave mixing and stimulated Raman scattering. We argue that an extremely strong 
nonlinearity of graphene in combination with its unique selection rules makes graphene 
a promising material for the new generation of compact optoelectronic devices. 



1.1. Band structure 

Graphene monolayer is a one-atom-thick monolayer of carbon atoms arranged in a 
hexagonal lattice, which we will treat as a perfect two-dimensional crystal structure 
in the (x, y)-plane. The electronic band structure of graphene has been extensively 
studied, starting from Wallace in 1947 [10] who used the tight-binding model. Here 
we briefly summarize the results relevant for our subsequent derivation of the optical 
response. The carbon atom in graphene has four valence electrons, three of which form 
tight bonds with the three neighbor atoms, therefore one atom only has one conduction 
electron in 2p z state. Considering only the interaction with the atom's three nearest 
neighbors, the resulting Hamiltonian in k-representation is purely off-diagonal: 
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Here 70 ~ 2.8 eV and a = 1.42 A are the nearest-neighbor hopping energy and C-C 
spacing. Then the energy dispersion relation is 

E k = ± lo f + 4cos ^ cos k -f + 4cos2 ^ (1) 

The electron and hole bands, denoted by ± are fully symmetric about the Dirac 
points at the six corners of the first Brillouin zone, 

v3k x a k y a 

— ^— = (2n + 1)tt, cos -|- = 0.5, 
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where = 0. Only two of these six Dirac points are inequivalent, referred to as K 
and K'. The Dirac points are located on the Fermi level if graphene is undoped and 
unbiased. In the vicinity of Dirac points, for example near k = K, it is convenient to 
define q = k — K . The effective Hamiltonian then becomes 

H 9 = hv F ( °. qx + iq y)=hv F ff-q, (2) 
\ q x -iq y J 

similarly to the one for an ultra-relativistic massless particle with spin 1/2 after replacing 
the velocity of light c with the band parameter (Fermi velocity) Vp = 370 /2ha ~ 10 6 
m/s. The pseudospin variable entering the problem is related to the presence of two 
sublattices A and B pQ. The corresponding eigenvalue is 



E@ = ±hv Fy /<£ + $, 
and the eigenfunctions are two-component spinors. 



1.2. Landau levels 

In an external magnetic field Bz perpendicular to the plane of graphene, the continuous 
energy bands near the Dirac points split into discrete Landau levels. The effective-mass 
Hamiltonian (TTJ H21 HBJ for a graphene monolayer can be written as a 4x4 matrix to 
combine the contribution of K and K' points: 
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where n = p + eA/c, p the electron momentum operator, and A is the vector potential, 
which is equal to (0, Bx) for a constant magnetic field. In this Hamiltonian the coupling 
between the K and K' point is neglected, so we can write down the solutions to the 
Schrodinger equation = e^fr separately for each point. For example, near the K 
point the Hamiltonian is Hk = Vpa-jf and the eigenfunction is specified by two quantum 
numbers n and k y , where n = 0, ±1, ±2, • • •, and k y is the electron wave vector along y 
direction: 
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where l c = \Jch/eB is magnetic length, H n {x) the Hermite polynomial. The eigen 

energy is 

e n = sgn(n)hu c ^J~\n\, u c = V2vp/l c . 

Positive or negative value of n corresponds to electrons or holes. Compared with 
Landau levels for a conventional 2D electron/hole system with a parabolic dispersion, 
E n = (n + 1/2) UeB/m*, Landau levels in graphene are unequally spaced: oc \/B. As 
shown in Fig. 1, the magnetic field "condenses" the original states in the Dirac cone 
into discrete energies, and each Landau level contains the same aerial density of states 
N§ = l/2irl c 2 , not including spin and valley degeneracy factors. 



2. Optical transitions between the Landau levels 

2.1. Selection rules 

Transitions between adjacent Landau levels in graphene fall into the mid-infrared 
to t erahertz (THz) range for a magnetic field in the range 0.01-10 Tesla: %oo c ~ 
36^B(Tesl a) meV. Consider an incident classical optical field E = E(u) exp (—iut)e 
polarized in the x-y plane along vector e. Let us define the left-hand circular polarization 
vector as clhs — [x — iy]/y/2 and the right-hand circular polarization vector crhs — 
[x + iy]/y/2. To include interaction with the optical field, we add its vector potential , 
A op t = icE/u, to the vector potential of the magnetic field in the generalized momentum 
operator n in the Hamiltonian. This results in adding the interaction Hamiltonian H int 
to ifo, where 

6 -* 

Hint = v F 3 ■ -A^ (5) 

Unlike the interaction Hamiltonian Hi nt for an electron with a parabolic dispersion, 
there are no higher order terms such as tt 2 near the Dirac point in graphene, so that 
even for a relatively strong optical field the interaction Hamiltonian is still linear with 
respect to A opt . Furthermore, H int does not contain the momentum operator; it is simply 
determined by the Pauli matrix vector a. The matrix element of the optical transition 
between Landau levels is given by 

(i\H int \j) = — (i\a x x + a y y\j) ■ E, 

where (i\a x x + cr y y\j) is 

V^CiCji-i)^^- 1 (sgn(rai) (0 K |_i|0| n .|) • e LHS + sgn(n J ) (^il^-i-i) • e RH s) ■ 

Since <p n are orthogonal, the above expression is nonzero only when \ni\ — 1 = \rij\ 
or |rij| = \nj \ — 1. As a result, the selection rule for the allowed transitions turns out to 
be 



(6) 
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where n is the energy quantum number. Denoting n/ and n, as the quantum numbers 
of the final and initial state, we can also conclude that crhs photons are absorbed when 
\rif \ = \rii\ — 1 while an absorption of a clhs photon leads to the transition |n/| = |r^| + l. 
Comparing with a typical selection rule for inter-Landau level transitions in a traditional 
2D system, An = ±1, the transitions with An greater than 1 are allowed in graphene, 
for example, from = —1 to n/ = 2, which leads to an efficient resonant nonlinear 
mixing. Mid/far- infrared optical absorption between Landau levels in monolayer and 
multilayer graphene has been extensively studied theoretically and in experiments; see 

e.g. |3l HI El Eld El 13- 

2.2. The dipole moment matrix of graphene 

To calculate the 2D optical polarization as an average dipole moment per unit area of 
the graphene sheet, 

P(F,t) = N(jt} = Ntr(p.fl, (7) 

where N is the surface density of electrons and p is their density matrix, we need to 
know the dipole moment matrix fx associated with inter-Landau level transitions. To 
calculate it, we first evaluate the commutator 

[r, H] = [r, Vpo ■ p\ + [r, vpcr ■ -A]. 

The second term on the right-hand side is zero because A is a function of f. So the 
commutator of f and the Hamiltonian is 

[r, H] = vpd • [f,p\ = ihvF<y. 

Since we use hats for both the operators and the unit vectors, we will stop putting hats 
over vector-valued operators unless it may create confusion. Choosing the eigen states 
of H as the basis, we obtain 

(m\[f,H]\n) = (m\fH\n) — (m\Hr\n) = (e n — e m )(rn\r\n), 

where e n and e m are the eigen energies of states |n) and |m). So the dipole matrix 
element of a closed system is defined as 

Pmn = e ■ (m\f\n) = — ^ — (m\v F a\n). (8) 

Similarly to the matrix elements of the interaction Hamiltonian (m\Hi nt \n) , the 
dipole matrix elements are determined by elements of the Pauli matrix (m\a\n). In 
particular, p mn is nonzero when \m\ = |n| ± 1. Using the wavefunction in Eq. (4), the 
analytic expression for the dipole moment element can be derived: 

tmn = ^rC* m C n (sgn(m)H)H-V| mhl , H)H0f m| ) 
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As an example, in a 4-level system that will be considered below (energy quantum 
numbers n= -1, 0, 1, 2), the allowed transitions are between n = ±1 and n = 0, and 
between n = ±1 and n = 2. The eigen functions of these four energy levels are 
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Combining with Eq. (|9|), the dipole moment matrix of the 4-level system is 

(11) 
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To summarize, the dipole moment of the transition between the Landau levels in 
graphene has a magnitude of the order of 



ehvp 



oc 1/y/B. 



This is a very large value for the transitions near the Dirac point where e n — e m ~ hu: c : 
vf/u c ~ 18 nm at B = 1 T. Note that the dipole moment grows rapidly, ~ A, with 
increasing transition wavelength. This is a faster growth than in atomic systems (~ \/A) 
or conventional semiconductors. Therefore, one expects a very strong nonlinear optical 
response in the mid-infrared and THz region. 

3. Linear optical response of graphene in a magnetic field 

3.1. Linear susceptibility of graphene 

The optical polarization of the graphene sheet, P(r, t) = Nti(p- /2), can be presented as 
a usual expansion in terms of electric susceptibilities if the density matrix is solved as a 
series in powers of the incident fields. For a weak monochromatic field it is enough to 
keep the first, linear in E, term in the expansion of the density matrix to find the linear 
susceptibility x '■ 
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The equation of motion for the density matrix elements is given by 

Pnm ~jT \E-n ^m) Pnm ~jT [Hi n f (t) , p\nm ^inm^yPnm Pnm ) " 

(12) 

Here we approximated incoherent scattering with phenomenological decay rates 7„ m 
describing the relaxation of the matrix elements to their equilibrium values p n e $; 
Pnm = for n 7^ m. 

Formal expansion of the density matrix in powers of the interaction Hamiltonian 
leads to the following differential equations: 

A(0) = _ • (0) _ / (0) _ leq)\ . 

Fnm ^nrnHnm i nm \rnm Hnm J i 

Pnm, = ~{^nm + lnm)P nm ~ jA_Hi n t, P ^]nm j 
Pnm = ~{} ijJ nm + lnm)P nm ~ ^flinti P ^}nm j 



where co nm = (e n — e m )/K. Choosing = p n e m \ we can calculate higher order terms 
step by step. The iteration formula is given by 

P { n2 = f ^[H m t{t'),p [N - l) ] nm ^V [fanm + Inm) • " t)]dt' . (13) 

The interaction Hamiltonian was derived above. It can be rewritten as 

e -» 

Hintit) = v F a ■ -A opt 

ievp _ -*. N 
= a ■ E(oj) exp (— mot) 

CO 

= - jl- E(u)exp(-iut) (14) 

Here we have defined 

~ ievF _ . ievp 

u = cr; (m m \n) = (m\a\n), 

to to 

which coincides with the dipole moment if the incident optical field is exactly on 
resonance with a given transition, that is p mn = pZ mn when e n — e m = Tluj. 

The first-order (linear) part of the density matrix can then be calculated from the 
iteration formula Eq. (13) and Eq. (14): 

Pnl = f ^[Hint{t'), P (eq) )n m enp [{lU nm + 7nm ) • (t' - t))dt' , (15) 

Jo n 

where 

[Hi n t(t ), ff- '']nm = (t^nvPvm ~ Pnv^ Pvni) ' E{t ) 

v 

= {Pmt-P^km-Eit')- (16) 

This yields the 2D first-order polarization in the form 

P w (uj)=NtT (p (1) /i) 



p (eq) _ p (eq) (fx ■ A 
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Here N is the 2D (sheet) electron density of graphene, which is n s n v N<s> = 2/(tc1 2 ) where 
n s = 2 and n v = 2 are spin and valley degeneracy For a left-hand polarized optical 
field, the circular polarization vector e is clhs = [x — iy]/\/2, and the term (U nm ■ ej 
in the above expression is nonzero only when |n| = |m| — 1. On the other hand, for a 
right-hand polarized optical field the term is nonzero only when |m| = |n| — 1. This of 
course corresponds to the polarization selection rules that were already derived above; 
see also [H El E]- Taking them into account, the expressions for the 2D linear optical 
susceptibility for the left /right-hand in-plane polarized optical field are: 



-AC 2 C 2 e 2 vl - />( c *) 
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3.2. Absorption coefficient 

The high-frequency absorbance in monolayer graphene at zero magnetic field, a = 
ire 2 /he, is a constant. In a high magnetic field, the absorption coefficient shows a 
series of peaks due to inter-Landau-level transitions. From the standard expression for 
a weak absorption, 

Mx (1) (w)], (19) 



and combining with Eq. (18), we can calculate the absorption coefficient of monolayer 



graphene for the left /right-hand in-plane polarized optical field: 

a(u,e LH s) = 
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The transition linewidth is of the order of 10 meV; for example, it was measured to be 
~ 30 meV in the magnetic field of 3 T [T3] . The corresponding relaxation rate 7 is then 
on the scale of a few ~ 10 13 s _1 . However, we should keep in mind that this number 
depends on the sample quality and the substrate used in the experiment. 

The above result agrees with the absorption coefficient calculated in [I] using the 
Keldysh's Green function approach. If we assume the relaxation rates between different 
levels to be the same, that is j nrn = 7, and follow their notation for Landau levels as 
(a, n), where n > and a = ±1 denote whether the corresponding state is in conduction 



(+) or valence (-) band, Eq. (20) can be rewritten in the form identical to the one in 

/a \ I'ihciaVn+luc-a'Jh^c) \P n > a ' Pn+l,a) 

a{u,e L HS) = 2^ ' 
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a(u,e RHS )= 




(21) 



., (a\/n — lu c — a'y/nuj c — u) 2 + 7 2 



4. Nonlinear optical response 

Strong optical nonlinearity of graphene, like most of its unique electrical and optical 
properties, stems from the peculiar energy dispersion of carriers near the Dirac points, 
E oc ±\p\. As a result, the electron velocity oc dE/dp induced by an incident 
electromagnetic wave is a strongly nonlinear function of induced electron momentum. 
Nonlinear electromagnetic response of classical charges with such an energy dispersion 
has been studied theoretically in [15]. Recently, the four- wave mixing in graphene 
without a magnetic field has been observed at near-infrared wavelengths [16]. Effective 
bulk third-order susceptibility was estimated to have a very large value, X® ~ io- 7 
esu, which is more than an order of magnitude larger than in gold films. 

Nonlinear cyclotron resonance in graphene was considered theoretically in [T7] . 
again in the classical limit, which can be applied only to electrons in a low magnetic 
field that occupy highly excited Landau levels n ^> 1, when energy and momentum 
quantization are neglected. In a recent work |5J we presented a quantum mechanical 
density-matrix description of the nonlinear optical response of graphene, which is valid 
for quantizing magnetic fields and strong optical fields, including the effect of saturation 
of inter-Landau level transitions. Due to unique optical selection rules for "massless" 
electrons near the Dirac point, one can implement a nonlinear interaction in which all 
optical fields are resonant to allowed optical transitions. The resulting magnitude of x^ 
turns out to be extremely large, of the order of 0.1 esu at mid/far-infrared wavelengths 
in the field of several Tesla. A similar strategy of a completely resonant nonlinear wave 
mixing has been implemented in asymmetric coupled quantum well systems, where 
one can increase the dipole moment of an intersubband transition involving a large 
change in the energy quantum number n by an appropriate band structure design 
[IEJ EH 1201 I2H 1221 123 However, the resulting third-order nonlinearity was still 
several orders of magnitude lower than in graphene for the same spectral range. 

4-1- Four-wave mixing 

Efficient nonlinear optical coupling becomes possible in graphene due to strong non- 
equidistance of the Landau levels, large magnitude of the dipole matrix elements, 
and unusual selection rules A|n| = ±1 which enable transitions with change in n 
greater than 1. In this section we study a specific example of the nonlinear optical 
interaction, namely the four-wave mixing. Consider a strong bichromatic field E = 
Ei exp(-iuit) + E 2 exp(—iu)2t) + c.c. normally incident on the graphene layer. Here U\ 
is nearly resonant with the transition from n = — 1 to n = 2 and E\ has left circular 
polarization. The frequency u 2 is nearly resonant with the transition from n = to 
n = ±1 and E 2 has linear polarization, so that it couples both to transition — 1 — > 
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and 1, as shown in Fig. 1. As a result of the four-wave mixing interaction, the 
right-circularly polarized field E 3 at frequency U3 = U\ — 2u2 nearly resonant with the 
transition from n = 2 to n = 1 is generated. 




Figure 1. Landau levels near the Dirac point superimposed on the electron dispersion 
without the magnetic field E = ±vp\p\. (b): A scheme of the four- wave mixing 
process in the four-level system of Landau levels with energy quantum numbers 
n = — 1, 0, +1, +2 that are renamed to states 1 through 4 for convenience. 

The frequencies involved in the four- wave mixing fall into the mid-infrared and THz 
region in the magnetic field of a few Tesla, as shown in Fig. 2. For example, at B — 3T, 
the nonlinear signal is generated at a wavelength of 48 fim in the presence of pump 
fields at wavelengths 8 and 20 fim. 




Magnetic Field (T) 

Figure 2. Transition frequencies in the 4-level graphene system shown in Fig. 1(b). 
LJ mn indicates the transition frequency between levels m and n. 



Truncating the master equation (12) to the 4-level system shown in Fig. 1(b) and 
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introducing slowly varying off-diagonal elements of the density matrix as pa = o^\e~ tU)lt : , 
P41 = c^e - ^ 3 *, and P32,2i = cr 32,2ie~ i< ^ 2 *,one can obtain the following set of equations 
for the amplitudes a nm in the steady state and in the rotating wave approximation: 



^21021 — ^21^21 — ^32^31 + ^41 cr 24 
«r 3 2(T32 = ^32^32 ~ ^43^42 + ^21°31 
^43(743 = ^4371,43 — ^41(713 + Vt\ 2 a 42 
«r 3 i(T3i = ^2lO"32 + ^4lO"34 ~ ^320"21 ~ ^43^41 
^42(742 = Vt*2iO-^i + ^32(^43 — VL/i2,a^2 ~ ^4lC"l2 

ir 41 cT4i = f2 41 ra 41 + f2 2 i042 — ^430-31. (22) 

Here the Rabi frequencies are defined as fiy = £?y and the population differences 

are = pa — p^-. The notation for the field amplitudes is as follows: £41 = Ei, E21 is 
the right circularly polarized component of E 2 , and £32 is the left circularly polarized 
component of E%. The complex dephasing T 4 i = 741 + i(co>4i —u>i) and similarly for other 
transitions; all detunings from resonance are small. 

If the incident field is not strong enough to perturb the populations, the population 



differences n mn are constant in Eqs. (22). As a result, the off-diagonal density matrix 
elements such as 043 can be solved analytically and written as an expansion in powers 
of the pump fields: 

^43 ^41^21^32 1741^21^32 ^41^21^32 

C43 = ^r^ n 43 - - 3r ™ ™ ™32 + . 3r ™ ™ n 2i + - 3r r r n41 

«r 43 ^r 43 r| 1 r| 2 ^i^r^r^ ^i^i^i^! 

+ J! 21 32 n 21 + 41 ' ^ n 4 i + - (23) 

The first term on the right-hand side describes the linear absorption and the next four 
terms describe the 3rd order nonlinear optical response; the higher-order terms are 
dropped. Note that the last term on the right-hand side corresponds to a stimulated 
Raman scattering of the pump field £41 into the signal field £43, which we consider in 
the next section. 

The optical polarization at the frequency U3 of the nonlinear signal in the rotating 
wave approximation is given by 



P(u 3 ) = N ■ a^e-^ 1 + c.c. 



Below we investigate different nonlinear effects contained in Eq. (23). Consider first the 



four- wave mixing interaction oj\ — 2ou 2 =^ ^3> described by the second through fourth 



terms on the right hand side of Eq. ( 23 ) 



Substituting the expression for 043 into P(w 3 ), and keeping only these three terms 
will lead to the third-order nonlinear susceptibility corresponding to the four-wave 
mixing: 

(3)/. . \ _ ^43^41/^32^21 ( P22 - P33 



(^) 3 r 43 V n^h 
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P22 ~ Pll _ Pll - P44 P22 - Pll \ , 

"P* "P* "P F V "T* I ^ ' 

1 31 1 21 1 42-L 41 J- 42-L 21 / 

To estimate the order of magnitude of x^ 3 \ we assume that all incident fields are 
in exact resonance so that the detuning factors = 7^ = 7 are real numbers and 
all dephasing rates are the same. We also assume for definiteness that state 1 is fully 
occupied while states 2, 3 and 4 are empty, which means pn = l,p 22 = P33 = P44 = 0. 
Then the expression for x^ is further simplified into 

(3)/ \ 3iV/i 43 /i4i/i32P21 

x M — ' (25) 

This expression contains a 2D electron density N and is a 2D (surface) susceptibility 
To convert it into the bulk susceptibility for comparison with other materials, we 
can divide it by the thickness of one monolayer Az ~ 3 A. Taking a reasonable 
value for the dephasing rate, 7 = 3 x 10 13 s -1 [H] , the bulk weak-field susceptibility 
xfl ~ 0.37(1/B(T)) esu = 5 x KT 9 (1/B(T)) m 2 /V 2 . Here the magnetic field is 
measured in Tesla. This is by far the strongest nonlinearity as compared to any material 
we know. 

When the incident fields increase in intensity, they start affecting populations on 



each level. In this case Eqs. (22) have to be solved together with the equations for 
diagonal components of the density matrix. Introducing phenomenological transition 
times Tij between levels i and j, we can write these equations as 

dm n 2 n 3 n A 

-TT = 7^ + 7^ + 7^ 1 ("21012 - "12021 + "41014 ~ "I404lj 

at 1 21 1 31 J 41 

dn 2 n 3 n 4 n 2 

—TT = 77, r — — 1 ("12021 ~ "21012 + "32 023 - "23032J 

dt 1 32 J 42 J 21 

dn 3 n 4 n 3 n 3 

—TT = 77 77 77 1 ("23 032 ~ "32023 + "43034 ~ "34043J 

at 1 43 i 3 2 J- 31 

n 1 + n 2 + n 3 + n 4 = 1. (26) 

This is of course a very crude approximation of the actual relaxation dynamics of 
electrons, but it allows us to estimate how the effects of the optical pumping and 
saturation affect the nonlinear mixing efficiency and power. 

It is convenient to normalize incident fields by their saturation values which 
determine the field strength at which the population at a given transition becomes 
significantly perturbed: 

hy/i2i/T 2 i n^ l32 /T 32 ftyW^i 

-^21 = 1 -^32= ! -^41 = " \ Z <) 

P21 p32 P41 

Then the corresponding saturation Rabi frequencies are given by Q nm = \Jj n m/T nm . 
We then introduce the dimensionless fields x, x', and y as 

x = n 21 /n s 2i; x' = n 32 /n s 32] y = n il /ni 1 . (28) 
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"^7 w x /(»,!/), (29) 



For estimation, we can simply assume that the relaxation rates T nm are the same, 
T nm ~ T. Then the solution to the density matrix equation of motion depends on the 
fields through only two dimensionless factors x and y. In particular, the scaling Eq. (25) 
for x i^s) becomes 

iV/i43/i4i/i32/i21 

W) 

where f(x,y) is a function of x and y shown in Fig. 3. It is equal to 3 when incident 
fields are weak, < 1, and quickly decreases as x and y become greater than one. 

The electric field of the generated signal is determined by the nonlinear polarization 
p( 3 \u3). From Maxwell's equations, neglecting the depletion of the pump fields, we can 
obtain 

dE 2ixu - , , 

= i P- 30 

oz c 

Note that here P is a 3D polarization (an average dipole moment per unit volume). 




Figure 3. Contour plot of f(x, y) as a function of normalized pump fields x and y. 



For a thin layer of graphene one can integrate Eq. (30) over the thickness of the 
layer and obtain 

M^) = i—X®WEi&) 2 , (31) 
c 

where is a 2D susceptibility. The magnitude of |^| grows with the pump at small 
pump intensities and decays at high intensities because of the decrease in x^- It reaches 
a maximum at x = 2.6, y = 1.56. Of course, these particular numbers depend on the 
relative values of the relaxation times between the Landau levels. However, the general 
conclusion that the maximum nonlinear signal is reached when the pump fields are of 
the order of the saturation values remains true. For fixed x, y ~ 1, x^ scales with 
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the magnetic field as B 1 , whereas E ~ E sat ~ \/i3. As a result, from Eq. (31), the 
maximum nonlinear signal scales with the magnetic field as 

\E™ X | ~ u 3X m M\E sat \ 3 ~ >/B±- • (^) 3 ~ R ( 32 ) 

£5 

If we define intensity as I = c\E\ 2 /87i, the intensity of the generated signal is related to 
the incident field intensities as 

hM = {^P) 2 \x {3) \ 2 hM(I 2 (u, 2 )f. (33) 

Fig. 4(a) shows the plot of I 3 as a function of the pump intensity I 2 when the 
second pump intensity ii is tied to I 2 by the optimal condition y = (1.56/2.6)x = 0.6a:. 
The conversion efficiency in the magnetic field of 1-10 T is I 3 /I 2 ~ 10~ 5 — 10~ 6 . This 




Figure 4. (a) Intensity of the 4-wave mixing signal as a function of the intensity of 
the pump field E 2 normalized by Iq — c|i? sa t| 2 /87r ~ 2.2 x 10 5 W/cm 2 . The value of 
Jo is the saturation intensity of the transition 1-2 calculated at B = 1 T and assuming 
that 1/T = 7 = 3 x 10 13 s _1 . I\ is set to satisfy y = 0.6a;. (b) Enlargement of (a) 
near the origin, which shows the intensity of the 4-wave mixing signal for a weak pump 
field. 

trend is reversed for a small incident pump intensity, when ~ 1/-B and f(x, y) is 
nearly a constant. As a result, for weak pumps -Z3 is higher in a smaller magnetic field, 
as illustrated in Fig. 4(b), which is the enlargement of Fig. 4(a) near the origin. 



4-2. Stimulated Raman Scattering 



The very last term on the right-hand side of Eq. (23) that we previously omitted 
describes another interesting nonlinear process: stimulated Raman Stokes scattering 
of the pump field Ei(u>i) into the field E 3 (oj 3 ); see Fig. 5. Note that this term does not 
depend on the second pump field ^2(^2); therefore for this section we can put E 2 = 0. In 
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this case the amplitude of the off-diagonal density matrix element 043, which determines 
the optical polarization at the frequency of the nonlinear signal, becomes 

I ^43 



043 



-^43 + 



|^4l| 2 ^43 



n 4 i / 1 + 1^41 17 (r 43 r; 



31, 



(34) 



v^W"* ■ i 3 r 43 r| 1 r| 1 

Here the complex detuning at the difference frequency is given by T 31 = 7 31 +i(u 3 i — uj\ + 
u 3 ), whereas other detunings are still T 41 = 7 41 + z(w 4 i — uj\) and T 43 = 7 4 3 + i(w 4 3 — co 3 ). 



4) 

31 



\2) 



1) 



Es 



n=2 



n=l 



n=0 



n=-l 



Figure 5. Raman Stokes scattering of the incident field E\ into the signal E3. 



Since the polarization P(u 3 ) is proportional to the field E 3 , the small-signal solution 
to the wave equation Eq. (30) has an exponential form, E 3 = E exp(gz), where g is 
given by 



2iru 3 N 3D fi 43 [i 43 



\n 



hcT 



n 4 3 



n\ 



13 



1 31 1 41 



n 4 i )/(l + 1^117(^3^)) 



(35) 



and N 3 d = 2/ (Aznl^) is the volume density of electrons in a layer of thickness Az. The 
real part of g gives the spectrum of the Raman gain. It is similar to the one derived for 
resonant Raman lasers in atomic and quantum- well systems [231 EH EH] . The gain peaks 
at the frequency of the two-photon resonance u>i — u 3 = lo 3 \. Its peak value increases 
when the pump frequency is tuned closer to the one-photon resonance U\ = u; 4 i. 

To estimate the maximum gain, we assume exact resonance for the pump and Stokes 
fields with corresponding transition frequencies w 4 i and a; 4 3, and take all dephasing rates 
to be the same, so that ~ 7. Then the gain factor is simplified to 



gAz 



4^3/4 



n 4 3 + 



\n 



r 



-n u 




(36) 



For a weak pumping |fi 4 i| 2 <C 7 2 , all population stays in the ground state of the system, 



1 and n i3 ~ 0. Then the maximum gain becomes 
gAz 



4^3/4, |^ 4 i| 2 



Ti^cll 7 2 



(37) 
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When expressed in usual dimensions cm/W for comparison with other materials, 
the Raman gain coefficient is really huge: around 20 cm/W in the magnetic field of 1 
Tesla, and assuming 7 = 3 x 10 13 s -1 . This is many orders of magnitude higher gain 
than the one reported for resonant intersubband Raman scattering in conventional 2D 
semiconductor systems: asymmetric coupled quantum well systems or quantum cascade 
lasers [23 EH ESI EZj ■ 




(a) X (b) 



Figure 6. (a) Population differences as functions of the normalized pump field, (b) 
Gain gAz for one monolayer of graphene as a function of the normalized pump field 
in the magnetic field B = 1 T. Solid line: the total gain, dashed line: only the Raman 
part assuming 71,43 = 0. 

For a stronger pump field, effects of the optical pumping and saturation become 



important. From the structure of the gain expression Eq. (36), it is clear that the gain 
reaches a maximum value when the pump field is of the order of the saturation value. 
This is a generic property of all resonant nonlinearities. For even higher fields, the gain 
drops due to a decrease in nu and an additional power-broadening term 1 + |f24i| 2 /7 2 
in the denominator. Using the same notation as in the previous section, we define the 
saturation field = hyj^jTu/ fai and the dimensionless pump field x = Ei/E^. 
Taking all relaxation times to be the same, ~ T, all population differences and the 
gain factor can be calculated analytically. They are shown in Fig. 6 as functions of the 
normalized pump field. Note that for our choice of equal relaxation rates, the optical 
pumping to the upper state 4 results in the population inversion on the signal transition: 
7143 > 0. This leads to an additional contribution to the gain, as is clear from comparing 
the total gain and the Raman contribution. The peak gain of about 2 % is amazingly 
high for just one monolayer of the material. By stacking several layers and placing the 
system in a high-Q THz laser cavity one can achieve a THz Raman laser with emission 
wavelength tunable by a magnetic field. 

In conclusion, we presented detailed studies of the linear and nonlinear optical 
response of graphene placed in a strong magnetic field. We showed that this system has 
an extremely high optical nonlinearity. We discussed two schemes of the nonlinear THz 
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generation in graphene based on the resonant four-wave mixing and Raman scattering of 
intense mid-infrared fields. The predicted nonlinear power makes graphene interesting 
for a variety of THz applications. Furthermore, one expects to find a similar physics of 
the nonlinear optical interactions in topological insulators, where the surface states have 
a massless dispersion and demonstrate a similar magneto-absorption on the transitions 
between the Landau levels [28J. 

This work has been supported in part by NSF Grants OISE-0968405 and EEC- 
0540832, and by the NHARP Project No. 003658-0010-2009. 
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